Diffuse interface approach to brittle fracture 
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We present a continuum model for the propagation of cracks and fractures in brittle materials. 
The components of the strain tensor e are the fundamental variables. The evolution equations 
are based on a free energy that reduces to that of linear elasticity for small e, and accounts for 
cracks through energy saturation at large values of e. We regularize the model by including terms 
dependent on gradients of e in the free energy. No additional fields are introduced, and then the 
whole dynamics is perfectly defined. We show that the model is able to reproduce basic facts in 
fracture physics, like the Griffith's dependence of the critical stress as a minus one half power of the 
crack length. In addition, regularization makes the results insensitive to the numerical mesh used, 
something not at all trivial in crack modeling. We present and example of the application of the 
model to predict the growth and curving of cracks in a non-trivial geometrical configuration. 

PACS numbers: 



I. INTRODUCTION 

There are many problems in condensed matter physics 
and materials science in which the aim is to describe 
sharp interfaces separating regions with qualitatively dif- 
ferent properties. This occurs for instance in solidifica- 
tion, dendritic growth, solid-solid transformations, grain 
growth, etc. Traditionally, two approaches have been 
followed to tackle these problems. In one of them, the 
problem is solved for a fixed position of the interface, 
and based on this, the expected evolution of the inter- 
face in the following time step is calculated, and the pro- 
cess repeated. This method has the practical drawback 
that the different structure of the interface at each time 
step makes necessary the full solution of a new problem 
each time. The second approach is sort of brute force, 
in which the system is modeled to the atomic scale, and 
evolved according to their (Newtonian) equations of mo- 
tion. The problem with this approach is that it is impos- 
sible in practice to span the many orders of magnitude 
between the atomic scale and the relevant macroscopic 
scale. 

The diffuse interface technique (including the so-called 
phase field models) is a novel powerful approach to study 
this kind of problems lj. It typically describes the sharp 
interface by an additional field <fr (or more than one) . For 
instance in a solidification problem <f> can be taken to be 
in the liquid and 1 in the solid. If the spatial variation of 
4> is known the interface can be located. Then the prob- 
lem of keep tracking of the interface is eliminated against 
having to include also as a new dynamical variable. <f> 
is coupled (usually phenomenologically) to the original 
degrees of freedom of the problem, and its dynamic evo- 
lution is not defined a priori, but has to be seeded from 
outside. 

A key ingredient in phase field models is regularization 
of the field <f>. Although the sharp interface is the sensible 
case, in order to implement the theory a smooth transi- 
tion of cj) between the values on both sides of the inter- 
face is necessary. Then the interface acquires a fictitious 



width, which however does not alter the physical behav- 
ior of the system if it is much smaller than any other rel- 
evant length scale. An additional, very important effect 
of regularization is to make the properties of the system 
independent of the underlying numerical mesh used to 
implement the problem on the computer. Regularization 
is usually achieved by the inclusion in the theory of terms 
that depend on gradients of <j), penalizing rapid spatial 
variations of this quantity. 

Within the field of fracture, the phase field models that 
have been proposed include those of Aranson, Kalatsky 
and Vinokur 0, Karma, Kessler and Levine 0, and 
Eastgate et al. (see also @)- All of them use an 
additional scalar field as a phase field, that is taken to 
be (asymptotically) zero within fractures and one within 
the intact material. 

There is now a general consensus that a complete de- 
scription of the fracture process cannot be given only in 
terms of macroscopic variables. In fact, the divergence 
of the stress field near to the crack tip implies that phys- 
ical conditions change a large amount on distances of 
the order of interatomic separation. Then details of the 
material at the atomic scale can have an effect on the 
macroscopic behavior of cracks. On the other hand, the 
roughly similar phenomenology of crack propagation ob- 
served in very different materials raises the expectation 
that a general description with a minimum amount of 
parameters dependent on microscopic characteristics is 
feasible. This is in the spirit of the phase field approach 
to fracture: one may think that the microscopic variables 
have an effect that translates in the form given to the en- 
ergy density of the phase field, in the form of the terms 
coupling the phase field to the elastic degrees of freedom, 
and in the dynamics assumed for it. Except in this ef- 
fective way, microscopic parameters do not appear in the 
phase field formalism. 

The phase field approach is already giving promising 
results. For instance, it has been shown that crack in- 
stabilities, oscillations and bifurcation can be obtained 
within this scheme @. The sharp interface limit of some 
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phase field models of stress induced instabilities has been 
studied in Its possible relevance to fracture is given 
in 0. 

We are going to present here a diffuse interface ap- 
proach that has some qualitative difference with previous 
ones. Most importantly, it does not introduce additional 
variables into the problem, the full set of variables are the 
components of the strain tensor £[9j. Description of frac- 
ture is achieved by the nonlinear form of the effective free 
energy density as a function of e. Actually, our energy 
is quadratic for small displacements (and then correctly 
describes linear elasticity) and saturates for very large 
displacements, then describing fracture (the saturation 
energy value being related to the fracture energy). Reg- 
ularization is provided by terms in the free energy of the 
generic form (de) 2 . 

There arc a number of reasons to pay attention to this 
model, both conceptual and from the point of view of 
implementation. First of all, the absence of additional 
degrees of freedom makes this model be probably the sim- 
plest continuous (non-atomistic) description of the frac- 
ture process. It is then interesting to know how, and to 
what extent, fracture phenomenology is captured by the 
model. 

From a practical perspective there are two important 
things to point out. First, an important characteristic is 
the tensorial nature of the variable describing the occur- 
rence of fractures. In the approaches in which a scalar 
field 4> is introduced, knowing that 4> has become zero at 
some point tells that a fracture is passing through that 
point, but does not tell anything about what direction 
the fracture follows. In our case fracture is described by 
e itself, and then if we know that a fracture is passing 
through some point, we can say immediately which di- 
rection it has. We believe this is an important computa- 
tional advantage: a single cell of a computational mesh 
is sufficient to encode the information about existence 
and direction of fracture. In models using a scalar field <p 
we need a whole neighborhood of a computational cell to 
encode the same information. Second, previous attempts 
to model fracture through non-linear elasticity have used 
typically the displacement field as fundamental variable. 
For those theories regularization is a problematic issue 
as it typically leads to higher order differential equations 
which are difficult to solve numerically. Our equations 
contain only second order derivatives of the strain field, 
and thus they much more smooth to solve numerically. 

In this paper we concentrate mainly in the presentation 
and validation of the model, giving only a short presen- 
tation of a non trivial application. We have divided the 
presentation in the following form. In the next Section 
we define the model and give analytically the structure 
on an infinite straight crack. Section III shows in detail 
that the model is able to reproduce the Griffith's crite- 
rion. In Section IV we present an example in which the 
crack paths are determined in a non trivial geometric 
configuration. In Section V we give a perspective of our 
planned future work with the model. 



II. THE BASIC MODEL 

The fundamental variables of the model are the com- 
ponents of the strain tensor e: 

where lij(r) are the displacements with respect to the 
unperturbed positions [l0| . Our approach follows closely 
that in Ref. ^1 ( see a l so E3) where it was used to 
study textures in ferroelastic materials and martensites. 
For clarity we present the model for a two dimensional 
geometry, the generalization to three dimensions being 
conceptually straightforward, although of course more in- 
volved. The symmetric tensor has three independent 
components. For convenience we will choose them in the 
form 

ei = (en + e 22 )/2 

e 2 = (e u -ea2)/2 (2) 

£3 = £12 = £21 

which are named respectively the dilation, deviatoric and 
shear components. These three variables are however not 
independent. In fact, since e is derived from the two 
displacements u\, U2, there is one constraint that has to 
be fulfilled, leaving only two independent variables. The 
constraint is known as the St. Venant condition and 
it can be written as 

{d 2 x + d>! - {d 2 x - d 2 y )e 2 - 2d x 8 v e 3 = 0. (3) 

It is easy to show that this is an identity if the definitions 
(JTJ and are used [3 ■ 

To define the model we need to know the form of the 
local free energy density Fl(e). To correctly describe 
elasticity for small displacements the limiting form of 
Fl for small e should be 



f l( £ ) = -jC l j k ie i je k i (4) 

where Ctjki is the fourth rank tensor of elastic constants 
of the material. We will specialize the expressions to an 
isotropic material, as this was our first aim to use this 
kind of model. In the isotropic case F^ reads 

Fl(e) = 2Be 2 1 + 2 f i(e 2 + e 2 ) (5) 

where B and [i are the two dimensional bulk and shear 
modulus of the material (related to the three dimensional 
values B 3D and /r 3D by fj, = /i 3D , B = B 3D + /i 3D /3 
for the case of plane strain, and fi = /i 3D , B = 
3B 3D ii/[B 3D + 4/i/3] for the case of plain stress). 
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The previous expression of the free energy must be ex- 
tended to large strains to account for fracture. To do this 
the energy must saturate for large deformation. This is 
the main requirement, and different choices can be made 
|15| . In the simulations presented below for isotropic ma- 
terials we have chosen the form 



1 + F° L (c)/fo 



(6) 



The limiting value fo of Fl for e — > oo is the energy 
density necessary to impose locally a very large value 
for at least one component of e, and it is then related 
to the crack energy of the problem (see below). From 
the present free energy form, we can say that a crack is 
nucleating when F® ~ fo, i.e., when typical values of e's 
are ej ~ {fo/B) 1 ^ 2 (assuming fi ~ B, as it is usually 
the case). Cracks in the system are thus characterized as 
regions where F® ^> fo. 

The second crucial ingredient of the model is regular- 
ization. That is provided by gradient terms in the free 
energy density. The terms we use are typically of the 
form 



a *( Ve *) 2 

i=l,2,3 



(7) 



with numerical coefficients aj. To retain rotational in- 
variance we have to choose a 2 = a 3 . In certain cases, 
and to avoid some unphysical behavior, these terms have 
to be cut off at large values of e» (see the next section 
for justification and details), and thus the gradient part 
of the free energy density is chosen as 



Fv = 



= 1,2,3 



aj(Vei) 2 s(ei) 



(8) 



where s(ej) goes to 1 for small ej, and tends to zero at 
large ej. The actual form of s(ej) we use is 



s{e l ) 



1 



i+(F° L /hr 



(9) 



where F® is given in JSJ, fx is a cut off value, and the 
exponent k controls the sharpness of the cut off (see next 
section). 

Once the full free energy is defined, the equations of 
motion are obtained by including a kinetic term T, which 
is quadratic in temporal derivatives of the displacements 
(i.e., T ~ 5/OU 2 , for some generic density p) and then 
it has to be transformed to a function of ej (see 0] 
for the details). The equations to be solved are more 
conveniently written in Fourier space (we use for the 
Fourier transforms of ej, and F = Fl + i*v). The result 
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where At, A 2l A3 are phenomenological damping coeffi- 
cients (that we will take to be equal Aj = A) and A is a 
Lagrange multiplier depending implicitly on E's, chosen 
to enforce at every moment the St. Venant constraint, 
that in Fourier space reads 



(*2 



kl)E 2 



2k x k y E 3 



0. 



11) 



In the examples below, the dynamic equations 1)10(1 are 
solved on a rectangular numerical mesh, using periodic 
boundary conditions in systems up to sizes of 512 x 512. 
In some cases we work with rectangular samples, with 
the length perpendicular to the crack being two or three 
times that along the crack, as we have observed that finite 
size effects are lower in this configuration than in a square 
sample with the same area. Also, the dynamic equations 
are solved in the overdamped regime, in which the second 
order time derivative terms in (|1(J|) are neglected. 



A. Infinitely long, straight crack. Analytical results 

As a necessary starting point, and also because it is 
probably the only case that can be treated analytically, 
we present here the structure of an infinite, straight crack 
in our model. In this geometry, all quantities depend 
only on the coordinate x (the crack is assumed to lie 
along the y direction), and the model becomes effec- 
tively one-dimensional. In fact, assuming the boundary 
conditions impose no strain in the y direction, we have 
et(x) = 62(2;), e3 = 0. As expected, there is a single 
order parameter in this configuration, that can be taken 
to be e\. The free energy of the model restricted to the 
present case becomes simply 



F/Ly = 



B*e 2 



dx + a |Vei| 2 s(ei)dx (12) 



B*e 2 /f 

where L y is the system size along the y direction, B* = 
2B + 2fi, a = a\ + a 2 , and the cut off function s is now 
a function of ei alone. For the choice in Eq. Q we get 
now 



s(ei) = 



1 



(13) 



\ + {B*e 2 /hY' 

It is clear that for a = the solution to Eq. i|12[l h as 
ex = except at a single point xq. The position xo is 
undetermined. This solution describes the system broken 
at xq . The fracture energy per unit length 7 in this case 
is non-zero only if a finite discretization 5 is used. In this 
case 7 = foS. 

To work out the finite regularization case (a / 0), it is 
more convenient to solve the problem under the assump- 
tion of an applied stress a along the x direction. To do 
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this we have to minimize the stress dependent free energy 
Fa, given by 



F a /L y = F/L y - 2cr J ei(x)dx 



(14) 



(the factor of two is due to the fact that in the present 
case e%= h^r)- Let us consider first the case fx — > oo, 
which implies (according to (|13fl ) s(ei) = 1, i.e., no cut 
off of the gradient terms at high values of ex- As in any 
one dimensional mechanical problem, the solution of Eqs. 
(|12H , (|14H is reduced to the evaluation of an integral. The 
numerical integration gives the profiles shown in Fig. ^ 
An important observation is that now the whole profile is 
smooth, and this implies that all macroscopic parameters 
that are numerically calculated will be independent of the 

mesh discretization S if this is small enough. 

For low er, in the region where ex ^>= \J fo / B* the 
following analytical solution is obtained 



ex(x) 



2a 



(15) 



where e™ ax is order fo/<r. The sharp fracture of the case 
a = transforms now into a smooth object that occupies 
a finite width w in the system. From (|15f) this width can 
be estimated to be 

w ~ V^/o 0- " 1 - (16) 

Note that w measures the width of the fracture in the 
original, unstrained reference system. An important pa- 
rameter to be calculated is the opening of the fracture A, 
defined as A = J (ex(x) — ei(oo)) dx. This has the mean- 
ing of the strain that the system is able to accommodate 
due to the existence of the crack. Its main contribution 
for low a comes from the central part of the crack, de- 
scribed by Eq. @5J. We get 



and from here and i|16fl 

,l/4 /o -l/4 Al /2 



a 



(17) 



(18) 



Then we see that upon increasing the opening A of the 
fracture, the width w increases as A 1 / 2 , and the stress 
decays only as A -1 / 2 . Then this model crack relieves 
the stress in the system only asymptotically. However, 
the decaying of the stress with A is too slow to give a 
finite crack energy. In fact, the crack energy per unit 
length 7 (defined as the energy of the system with the 
crack, minus the energy of the system at the same stress 
without the crack) can be estimated as: 
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f 1/4^3/4.1/2 

f w ~ a 1 f ' A ' 



(19) 



This means that the present regularization scheme does 
not provide a crack with a finite energy in the limit of 
large stretching, i.e., A — > oo. 



In the next section we will argue that this behavior 
does not invalidate completely the use of the present al- 
gorithm with fx — > oo if we are interested in cracks that 
do not become infinitely long. Notwithstanding, if we 
want a more accurate description of cracks, an implemen- 
tation in which the energy density of the crack remains 
finite as its length diverges is mandatory. The following 
is a possibility to obtain this. 

In order to have a model that relieves more efficiently 
the elastic energy, we first remind that gradient terms are 
necessary for a fracture to propagate without interference 
of the numerical mesh, therefore, they are important only 
when fracture is forming, namely, when typical values of 
ei are close or lower than ~ (fo/B) 1 / 2 . In the regions 
in which the fracture has nucleated, typical values of e's 
become much larger, and the gradient terms are not nec- 
essary any more. Then we can weaken their effect in 
those regions by choosing a non-trivial cut off function 
as given by Eq. © (or Eq. (fH5f in the present one- 
dimensional case). It should be emphasized that this ad 
hoc modification of the free energy does not introduces 
new parameters in the relevant regions that determine 
crack growth, namely, close to the crack tips. 

Our aim is to choose a value of n that generates cracks 
with finite width and energy in the limit A — » oo. I n or- 
der to do this, we first notice that for ex ^> \J fo / B* (i.e., 
where elastic forces become vanishingly small), the Eulcr- 
Lagrange equation corresponding to Eqs. (|12fl . (|13|) . and 
(O leads to 



1 + 



(¥ 



de\ 
dx 



— — ) + 2oe\ = constant 



(20) 



In particular, if fx — ► oo, we re-obtain from here the 
behavior given in fL5l) . The constant in Eq. ^ 
must b e calcu lated matching the pr esent solution for 
e i 3> yfo/B*, with that for ex < y/fo/B*. Assuming 
/i 3* fa for simplicity, it is obtained that this constant 
is Dfo, where D is a (order 1) numerical factor. Upon 
integration, Equation l|2U|) allows to find the full profile 
of the crack (in Fig. [21 we can see the results of numer- 
ical integration for k — 2). We concentrate in the case 
a — ► 0, in which the crack has relaxed completely the 
applied stress. The crucial result is that in this case, and 
for K > 1, it exists a well defined profile of the crack, 
given by 



(21) 




where the function g is defined as 

dw 



9W 



In VI + w 2k 
For x — > the limiting form is 



(22) 



(23) 
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The opening of the crack A = J e±(x)dx gives a finite 
value if k > 2. This means that the system has com- 
pletely relaxed the applied stress with a finite opening 
of the crack. Any further increase A* of A is accommo- 
dated in the system at the center of the crack, through 
a singular term A*S(x) added to f2H|l . For 1 < k < 2, 
the value of A is divergent due to the non-integrable di- 
vergence of ei around the origin in i|23[l . In this case 
the stress is relaxed only asymptotically, but still rapidly 
enough (compared to the /i — + oo case) to guarantee a 
finite width w and energy 7 of the fracture. In fact, from 
(|23[) the width w of that part of the system for which 
e i ~ yfo/B* can be estimated to be 

in the same way, the energy of the crack 7 becomes finite, 
and its value is proportional to 

, If aha , . 

7 ~ fow ~ y (25) 

From a numerical point of view, the algorithm with k > 
2 is probably too singular to be implemented successfully. 
We have implemented the case k — 1.5, which we have 
found is numerically tractable, and provides and almost 
perfect verification of the Griffith's criterion even in the 
quite small system sizes that we are using. 



III. FINITE STRAIGHT CRACK UNDER 
MODEL I LOADING: VERIFICATION OF THE 
GRIFFITH'S LAW 

One of the basic cornerstones of fracture physics, de- 
scribed in the very first pages of any fracture book, is the 
so called Griffith's criterion (GC). It states that under 
the application of a remote stress cr on a (effectively two- 
dimensional) system with a preexistent crack of length I, 
the crack will extend and eventually break the sample if 
cr is larger than some critical value a cr which scales as 
l~ x / 2 . A simple justification of this behavior can be given 
on energetic grounds. Upon the application of the stress 
cr, the system with the crack stores an elastic energy that 
is reduced in a quantity E e i = Ca 2 l 2 with respect to that 
of the system without the crack (C is proportional to an 
elastic constant of the material). On the other hand, 
the creation of a crack of length I is assumed to have an 
energy cost of E cr — jl, where 7 defines the fracture en- 
ergy per unit length. The total energy of the system as 
a function of I is then given by 

E = -E el + E cr = -Ca 2 l 2 + jl (26) 

and it has a maximum as a function of I at l m ax = 
7/(2Ccr 2 ). If I > l max , the system relieves sufficient 
elastic energy upon crack length increase to pay for the 



energy cost of crack creation, and the crack typically ex- 
tends abruptly and breaks the material. If I < lmax there 
is no sufficient energy in the system to make the crack 
enlarge. Note that according to H26[). the system would 
prefer to reduce I in order to reduce its energy. Experi- 
mentally this 'healing' of the crack is prevented by irre- 
versible processes: the crack energy is not recovered upon 
crack healing, and then the situation is that any I < l max 
is a stable crack. The previous arguments do not depend 
on the orientation of the crack in the system, assuming 
the material and the remote applied stress are isotropic. 
The GC is at the base of the fracture of brittle materials, 
and has to be satisfied by any model devised to describe 
such process. 

In our simulations we control the mean uniform strain 
e\. This correspond to an isotropic stress cr = 2Be\. In 
order to correctly calculate the critical stress cr c for cracks 
of a given length, we start from an initial spatial distri- 
bution of the variables e±, e2, e% that roughly describes a 
crack, and apply a stabilization procedure through a neg- 
ative feedback loop in the program, that monitors the 
length of the crack (defined using the contour defined 
by the value of the elastic energy F®=2B), and reduces 
(increases) the applied strain when the length increases 
(decreases). Results in Figs. |21 El and [7| below, where 
obtained with this procedure. 

The model with no regularization satisfies the GC if 
we restrict to cracks running in a single direction with 
respect to the underlying numerical mesh. This is shown 
in Fig. |3J We also see that no noticeable system size 
effects are observed for the system sizes used. As we use 
periodic boundary conditions this means that the crack 
is not influenced by the elastic field of its neighbor im- 
ages. However, the value of a c is strongly dependent on 
the orientation of the crack as shown in Fig0] This is 
one typical drawback of many discrete models of fracture 
when trying to simulate isotropic materials. Moreover, 
the influence of the numerical mesh has a clearly visible 
manifestation: as Fig. shows, if a crack is placed at a 
finite angle with respect to the mesh, when the critical 
stress is reached, the crack opens along one of the lattice 
main directions, instead of extending along its original 
direction as it should do in an isotropic system. 

The gradient terms are included in the model to solve 
this problem, and to make the behavior isotropic. The 
question of the satisfiability of the GC in the presence 
of a non-zero oti is not trivial since as we showed in the 
previous section, the energy of the crack per unit length 
does not saturate upon increasing strain unless an ade- 
quate cut off of the gradient terms is included. Let us first 
study the effect of a finite regularization without cut off 
(fi — ► 00 in Eq. ©). The fact that in this case the crack 
energy per unit length of an infinite crack is divergent (as 
seen in the previous section) is a manifestation of the fact 
that the crack energy of a crack of finite length I grows 
more rapidly than I itself. On the basis of the energetic 
arguments for the GC, we expect here a dependence of 
cr c on I of the form cr c ~ with [3 < 1/2. 
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Numerical simulations first of all confirm that the be- 
havior of the system can be made isotropic by including 
regularization. In fact Fig. Q]shows that the critical stress 
become independent of the angle between the crack and 
the numerical mesh. More than that, the crack extends 
along the original direction it had, independently of the 
numerical mesh (see Fig. EJ. However, as anticipated, 
a power low decaying with an exponent lower than 1/2 
is obtained for the critical stress as a function of crack 
length (see Fig. 0). The fitted power for the parameters 
used is ft = 0.402 =p 0.009. Then the model with finite 
on but infinite /i gives a slightly incorrect behavior of 
critical stress as a function of crack length. If this dis- 
crepancy can be considered small, then the model with 
infinite fx (which is easier to implement) can be used. If 
the previous discrepancy is considered serious (whether 
it is serious or not will depend on the particular problem 
under study) then a finite f\ formalism has to be imple- 
mented. As already discussed, a power k > 1 in Eq. © 
has to be used to guarantee that the GC is satisfied. In 
Fig. 0we show results indicating that a very good fitting 
to the GC is obtained for n — 1.5 and /i = 14.65. This 
value of /i was chosen to obtain the best verification of 
GC in the finite system we are using, but in principle 
any (finite) value of /i should reproduce the 1/2 power 
dependence in the case of sufficiently large system sizes. 
Fig. [8] shows the three dimensional profiles of the crack 
with I = 83(5, stabilized at the critical stress. 



IV. ELASTIC INTERACTIONS BETWEEN 
CRACKS 

The prediction of crack trajectories under general load- 
ing conditions for bodies of arbitrary shape and possibly 
with pre-existent cracks is very important in many engi- 
neering applications. The present formalism is well suited 
to study this kind of problems, in particular in those cases 
in which slight deviations from straight propagation are 
expected. These cases are particularly difficult to tackle 
with a non-regularized model. As a very simple an il- 
lustrative example of that, we consider a pair of parallel 
cracks loaded isotropically. We show here some qualita- 
tive results, leaving a more detailed quantitative analy- 
sis for a forthcoming work. We first stabilized the two 
parallel cracks by the feedback mechanism already ex- 
plained, then stop this stabilization, and increase a small 
amount the stress, and follow the crack evolutions as they 
grow. Snapshots of the system during crack propagation 
(Figs. 1 1 U|) show that the cracks propagate diverging 
from straight propagation. This is a non-trivial effect 
caused by the elastic interaction between the two cracks. 
Note that in the present case, cracks propagate from both 
fractures, because of the perfect mirror symmetry of the 
problem with respect to the middle plane. 

In a slightly different configuration, we place the two 
parallel cracks shifted (Fig. 111(1 . Now the propagation 
of cracks from the internal tips is strongly influenced 



by the nearby crack, producing a geometrical pattern of 
curved cracks that is well known (see for instance (19j). 
The propagation from the external tips is now almost 
not influenced by the second crack and then essentially 
straight. 



V. SUMMARY, OUTLOOK AND 
CONCLUSIONS. 

In summary, we have presented a model for the study 
of cracks propagating in brittle materials. The model 
does not use additional variables others than the strain 
tensor, and then is a minimalist continuum model for 
description of cracks. As a crucial ingredient it includes 
regularization terms that make the cracks be smoothed. 
We have shown that in this form the model can describe 
accurately an isotropic material in conditions in which 
the non-regularized model fails neatly. 

We have validated the model showing how it can fit 
accurately the Griffith's law for the critical stress as a 
function of crack length. We have seen that in order to 
obtain this result the regularization has to be softened 
in the interior of the cracks. This ad hoc modification 
however leaves intact the model in the neighborhoods of 
the crack tips, where the processes responsible for crack 
advance take place. As a first application we have shown 
how the model can predict the propagation and eventual 
bending of cracks induced by elastic interactions between 
them. 

We believe that the present technique is straightfor- 
ward to implement and computationally efficient, and 
that it addresses in a phenomenological way the very im- 
portant applied problem of predicting crack evolution, 
without requiring as explicit input any details about the 
physical conditions in the process zone. The technique 
can be implemented also in three dimensions, although 
we feel that this should wait for some increase in com- 
putational power before this can be implemented on a 
desktop computer. 

We want to indicate a few important direction along 
which the model can be applied, and in which we have 
started some work. First of all, all simulations in the 
present paper were done in the overdamped regime, 
where dynamical effects are absent. This may be a rea- 
sonable choice for the cases in which the cracks are known 
to grow quasistatically. This may include crack propaga- 
tion under slowly varying external conditions, as for in- 
stance non-uniform thermal stresses. For cases in which 
dynamical effect are expected to be important the full dy- 
namical equations have to be implemented. Preliminary 
work shows that indeed the implementation of the iner- 
tial dynamics gives rise (under appropriate conditions) 
to well known phenomena such as crack oscillation and 
bifurcation. We expect to report about this soon. An- 
other possible interesting application of the model con- 
cerns the determination of minimum energy configura- 
tion of cracks. In fact, as a result of regularization, our 
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cracks are in principle able to shift laterally, in addition 
to extend from its tips. This effect is not observed in 
the simulations presented here as it occurs typically in 
much longer time scales than the one we were interested 
in, but it can be enhanced under particular conditions. 
This shifting of cracks is driven by the tendency of the 
system to minimize its energy, and then it provides a tool 
to study cases in which the minimum energy configura- 
tion has a physical meaning. An example of this kind of 
application has been presented in [2IJ 

As a final consideration, we stress that we are describ- 
ing fracture in a continuum model of a brittle material 
as a non-linear elastic process. Due to the simplicity of 
the model, the influence of microscopic details at the pro- 



cess zone have only very few places were to leak in the 
present formalism. One point where this may happen 
is in the form of the interpolation function between lin- 
ear elasticity and broken material regime. Based on very 
recent findings. |20l| we expect that particular changes in 
the form of this function may give rise to different phe- 
nomenological behavior of crack propagation. 
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FIG. 1: Profiles of ei for the one- dimensional problem, at dif- 
ferent values of the global strain ei = SL X /L X . Labels are the 
values of (B*//o) 1/ ' 2 ei used for each curve. Note that ei does 
not go to zero away from the fracture, but a remmanent strain 
remains. This value however, decreases when the fracture is 
opened wider. 
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FIG. 2: Same as previous figure when the gradient terms 
are cut off at large values of ei, as explained in the text. 
Parameters used are /i//o = 100, k = 2. Note that contrary 
to the case in the previous figure, the fracture width (and this 
implies also the fracture energy) saturates to a finite value for 
large ei. 
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FIG. 3: (Color online) Critical remote stress for the propa- 
gation of straight cracks of different lengths I (with 5 being 
the lattice discretization), placed along the x direction, for 
two different system sizes (a value of /i = 0.5B is used in all 
numerical simulations presented). Results are for the model 
without regularization (a — 0). The exponents to fit the 
Griffith's law (a c ~ with /3 = 1/2) are indicated. 
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FIG. 4: (Color online) Critical remote stress a c for differ- 
ent values of a' = a/2B5 2 , normalized by the mean values, 
{(T C ) = 0.2436B (for a' = 0), (<r c ) = 0.9254B (a' = 0.25), 
(<r c ) = 1.0776.B (a' = 0.5), for the propagation of cracks of 
length la = 405 placed at different angles with respect to the 
numerical mesh (256 x 256). The strong dependence with an- 
gle in the a — becomes barely noticeable after the inclusion 
of regularization. 
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FIG. 5: (Color online) Effect of the numerical mesh on the 
propagation of cracks in the model without regularization 
(a = 0). Snapshots of cracks of original length lo = 165, 
growing with 9 = Odeg (upper plots), and 6 = 30deg (lower 
plots) and a remote isotropic strain ei = 0.19, larger than 
the corresponding critical value. The left panels are the con- 
figuration soon after the application of the loading, and the 
right panels are those some time later. The influence of the 
numerical mesh is evident. The full simulated system size is 
256 x 256. The key to the gray (green) scale here and in all 
following figures is as follows: We plot the values of ei from 
brighter (lowest ei) to darker color for a value ei>2 (cor- 
responding to the crossover to cracked material). All values 
above that one are plotted black and correspond to the region 
inside the crack. 




FIG. 6: (Color online) Same as previous figure for a = 0.5BS 2 , 
/i — ► oo, lo = 408 and ei =0.56 . The effect of the numerical 
mesh is undetectable. 
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FIG. 7: (Color online) Critical remote stress for the propaga- 
tion of straight cracks of different lengths I, obtained including 
into the model the regularization terms, a' = 0.5. An expo- 
nent lower than 0.5 is obtained with infinite /i but a /3 = 0.5 
exponent is well fitted with a finite /i value. System size: 
192 x 512. 



0-3 r 




FIG. 8: The full profiles for the three variable ei, e2, and 
eg of a crack with / = 835, at the critical stress a c = 0.7B, 
corresponding to the encircled point in Fig. [7| 
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FIG. 9: (Color online) Time sequence for two parallel cracks 
(length lo — 60S, separation equal to 205) loaded isotropically 
under ei= 0.516 (larger than the corresponding critical value), 
with a = BS 2 , fi — 14.65, K = 1.5. The cracks extend 
with some divergence, due to a non-trivial effect of elastic 
interaction between them. 



FIG. 10: Time sequence in contours of the elastic energy at 
F° = 2B for the parallel cracks of the previous figure (equally 
spaced each At = 25A/B, arrows indicate increasing time) 
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FIG. 11: Same as previous figure for skew-parallel cracks 
with initial length lo = 100(5, perpendicular distance between 
them equal to 50<5 and horizontal distance between the outer 
tips equal to 2005 (simulated system size: 384 x 512). The 
isotropic load is ei=0.3, corresponding to a stress a = 0.6B, 
larger than the corresponding critical value a c ~ 0.5B. 



